/*
 * CompareDNA.cpp
 *
 *  Created on: May 3, 2011
 *      Author: HTPC
 */
#include <iostream>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
//#include <windows.h>
#include "DNA.h"
#include "CompareDNA.h"
#include "FileIO.h"

using namespace std;

//Initializing the member variables
vector< vector< int > > ComparisonMatrix;
vector< vector< int > > ScoringMatrix;
vector< vector< int > > DistanceMatrix;
vector< vector< int > > AdjustedDistanceMatrix;

int openingGapPenalty;
int extendingGapPenalty;
double Matches;
int percentDifference;

CompareDNA::CompareDNA(){

}

void CompareDNA::setValueInDistanceMatrix(int i, int j){
	DistanceMatrix[i][j] = percentDifference;
}

void CompareDNA::initializeDistanceMatrix(int size){
	DistanceMatrix.resize(size);
		for(int i = 0; i < size; i++){
			DistanceMatrix[i].resize(size);
		}
}
///Creates the Newick UPGMA
void CompareDNA::buildNewickUPGMA(vector<string> FileList){
	joiningUPGMA(FileList);
}

///Is responsable for the creation of the newick tree text
void CompareDNA::joiningUPGMA(vector<string>FileList){

	string NewickString = "";
	double leftLegWeight = 0;
	double rightLegWeight = 0;
	double legLength;
	double leftLeg;
	double rightLeg;

	NewickString = '(';
	//begin loop
	vector<int> SmallestNumber = findSmallestDistance();
	legLength = DistanceMatrix[SmallestNumber[0]][SmallestNumber[1]] / 2;

	leftLeg = legLength-leftLegWeight;
	rightLeg = legLength-rightLegWeight;

	char Integer[33];
	itoa(leftLeg, Integer, 10 );
	NewickString += '(' + FileList[SmallestNumber[0]] + ':' + Integer + ',';
	itoa(rightLeg, Integer, 10 );
	NewickString += FileList[SmallestNumber[1]] + ':' + Integer + ')';

	leftLegWeight += leftLeg;
	rightLegWeight += rightLeg;

	//initialize size of the adjustedDistanceMatrix
		AdjustedDistanceMatrix.resize(DistanceMatrix.size()-1);
		for(int i = 0; i < (int)DistanceMatrix.size()-1; i++){
			AdjustedDistanceMatrix[i].resize(DistanceMatrix.size()-1);
		}


	//endloop

	NewickString += ')';

	//update leftlegweigth and rightlegweight
}

///Finds the Distance in a Distance Matrix
vector<int> CompareDNA::findSmallestDistance(){
	int smallestDistance = 0;
	vector<int> coordinates(2);

	for(int i = 0; i < (int)DistanceMatrix.size(); i++){
		for( int j = 0; j < (int)DistanceMatrix[0].size(); i++){
			if(j > i && j != i){
				if(smallestDistance > (int)DistanceMatrix[i][j]){
					smallestDistance = (int)DistanceMatrix[i][j];
					coordinates[0] = i;
					coordinates[1] = j;
				}
			}
		}
	}
	return coordinates;
}

///Compares to pieces of DNA and returns there score for making a scoring matrix
int CompareDNA::CompareResidue(char DNAResidue1, char DNAResidue2){

	int i = 0;
	int j = 0;

	if(DNAResidue2 == 'A')
		i = 0;
	if(DNAResidue2 == 'G')
		i = 1;
	if(DNAResidue2 == 'C')
		i = 2;
	if(DNAResidue2 == 'T')
		i = 3;

	if(DNAResidue1 == 'A')
		j = 0;
	if(DNAResidue1 == 'G')
		j = 1;
	if(DNAResidue1 == 'C')
		j = 2;
	if(DNAResidue1 == 'T')
		j = 3;
	return ScoringMatrix[i][j];
}

///Calculates the percent difference in to DNA sequences
void CompareDNA::calcPercentDifference(DNA *DNASequence1, DNA *DNASequence2){
	//Initialize Scoring Variables
	MisMatch = 0;
	percentDifference = 0;
	for(int i = 0; i < DNASequence2->getAlignedSize()+1; i++){
		if(DNASequence1->getDNAatPositionAligned(i) != DNASequence2->getDNAatPositionAligned(i)){
			MisMatch++;
		}
	}
	percentDifference = (MisMatch/DNASequence2->getAlignedSize()) * 100;
	cout << "Percent Difference: " << percentDifference << endl;
}

///Prepares the comparison matrix to be used
void CompareDNA::initializeComparisonMatrix(DNA *DNASequence1, DNA *DNASequence2, int openGapP){

	/*If a sequence has already been aligned using a different method, the
	 * aligned sequence needs to be initialized as well.
	 */
	DNASequence1->clearAligned();
	DNASequence2->clearAligned();
	//Clearing the matrix just in case something else was there
	ComparisonMatrix.clear();
	//setting the size of the two dimensional vector
	ComparisonMatrix.resize(DNASequence2->getOriginalSize()+1);
	for(int i = 0; i < DNASequence2->getOriginalSize()+1; i++){
		ComparisonMatrix[i].resize(DNASequence1->getOriginalSize()+1);
	}

	//intialization with gap penalty
	openingGapPenalty = openGapP;
	ComparisonMatrix[0][0] = 0;
	for(int i = 1; i < DNASequence1->getOriginalSize()+1; i++){
		ComparisonMatrix[0][i] = ComparisonMatrix[0][i-1] + openingGapPenalty;
	}
	for(int i = 1; i < DNASequence2->getOriginalSize()+1; i++){
		ComparisonMatrix[i][0] = ComparisonMatrix[i-1][0] + openingGapPenalty;
	}
}

///Initializers for the scoring matrix
void CompareDNA::basicScoringMatrix(){

	ScoringMatrix.clear();
	ScoringMatrix.resize(4);

	for(int i = 0; i < 5; i++){
		ScoringMatrix[i].resize(4);
	}

	ScoringMatrix[0][0] = 1;
	ScoringMatrix[0][1] = -1;
	ScoringMatrix[0][2] = -1;
	ScoringMatrix[0][3] = -1;

	ScoringMatrix[1][0] = -1;
	ScoringMatrix[1][1] = 1;
	ScoringMatrix[1][2] = -1;
	ScoringMatrix[1][3] = -1;

	ScoringMatrix[2][0] = -1;
	ScoringMatrix[2][1] = -1;
	ScoringMatrix[2][2] = 1;
	ScoringMatrix[2][3] = -1;

	ScoringMatrix[3][0] = -1;
	ScoringMatrix[3][1] = -1;
	ScoringMatrix[3][2] = -1;
	ScoringMatrix[3][3] = 1;
}

///Does a standard scoring method
void CompareDNA::standardScoringMatrix(){
	ScoringMatrix.clear();
	ScoringMatrix.resize(4);
	for(int i = 0; i < 5; i++){
		ScoringMatrix[i].resize(4);
	}

	ScoringMatrix[0][0] = 10;
	ScoringMatrix[0][1] = -1;
	ScoringMatrix[0][2] = -3;
	ScoringMatrix[0][3] = -4;

	ScoringMatrix[1][0] = -1;
	ScoringMatrix[1][1] = 7;
	ScoringMatrix[1][2] = -5;
	ScoringMatrix[1][3] = -3;

	ScoringMatrix[2][0] = -3;
	ScoringMatrix[2][1] = -5;
	ScoringMatrix[2][2] = 9;
	ScoringMatrix[2][3] = 0;

	ScoringMatrix[3][0] = -4;
	ScoringMatrix[3][1] = -3;
	ScoringMatrix[3][2] = 0;
	ScoringMatrix[3][3] = 8;
}
//End Initializers for the scoring matrix

///Prints the scoring matrix to the screen
void CompareDNA::printScoringMatrix(){
	for(int i = 0; i < 4; i++){
			for(int j = 0; j < 4; j++){
				cout << ScoringMatrix[i][j] << " ";
			}
			cout << endl;
		}
}

///Prints the comparison matrix to the screen and writes it to the file
void CompareDNA::printComparisonMatrix(){
	//printing comparison matrix
	for(int i = 0; i < (int)ComparisonMatrix.size(); i++){
		for(int j = 0; j < (int)ComparisonMatrix[0].size(); j++){
			printf(" %*d ", 5,ComparisonMatrix[i][j]);
            FileIO::Write(ComparisonMatrix[i][j]);
		}
		cout << endl;
        FileIO::Write('\n');
	}
}

///Prints the distance matrix to the screen and writes it to the file
void CompareDNA::printDistanceMatrix(){
	for(int i = 0; i < (int)DistanceMatrix.size(); i++){
			for(int j = 0; j < (int)DistanceMatrix[0].size(); j++){
				printf(" %*d ", 5,DistanceMatrix[i][j]);
	            FileIO::Write(DistanceMatrix[i][j]);
			}
			cout << endl;
	        FileIO::Write('\n');
		}
}

///Prints two aligned sequences to the screen and writes them to a file
void CompareDNA::printTwoAlignedSequences(DNA *Sequence1, DNA * Sequence2, string type){
	cout << "Enter the width of characters you would like to print on each line:" << endl;
	int width;
	try{
		cin >> width;
	}catch(exception){
		cout << "An integer is needed...try again...exiting" << endl;
	}

	cout << "Two DNA Sequences that have been aligned using: " << type << endl;
	FileIO::Write("Two DNA Sequences that have been aligned using: ");
	FileIO::Write(type);

	//HANDLE hConsoleOut; //handle to the console
	//hConsoleOut = GetStdHandle(STD_OUTPUT_HANDLE);
	int DNA1Position = 0;
	int DNA2Position = 0;
	int i = width;
	cout << "printing out DNA sequences line by line for easy comparison:" << endl;
	FileIO::Write("printing out DNA sequences line by line for easy comparison:\n");
	cout << endl;
	do{
		/*while(DNA1Position <= i && DNA1Position < Sequence1->getAlignedSize()){
			if(Sequence1->getDNAatPositionAligned(DNA1Position) == '_'){
				SetConsoleTextAttribute(hConsoleOut, FOREGROUND_RED);
			}else{
				SetConsoleTextAttribute(hConsoleOut, FOREGROUND_BLUE | FOREGROUND_RED | FOREGROUND_GREEN);
			}*/

			printf("%c", Sequence1->getDNAatPositionAligned(DNA1Position));
			FileIO::Write(Sequence1->getDNAatPositionAligned(DNA1Position));
			DNA1Position++;
		//}
		cout << endl;
		FileIO::Write("");
		/*while(DNA2Position <= i && DNA2Position < Sequence2->getAlignedSize()){
			if(Sequence2->getDNAatPositionAligned(DNA2Position) == '_'){
				SetConsoleTextAttribute(hConsoleOut, FOREGROUND_RED);
			}else{
				SetConsoleTextAttribute(hConsoleOut, FOREGROUND_BLUE | FOREGROUND_RED | FOREGROUND_GREEN);
			}*/
			printf("%c", Sequence2->getDNAatPositionAligned(DNA2Position));
			FileIO::Write(Sequence2->getDNAatPositionAligned(DNA2Position));
			DNA2Position++;
		//}
		i = i + width + 1;
		cout << endl;
		cout << endl;
		FileIO::Write("\n");
	}while(DNA1Position < Sequence1->getAlignedSize());

	//SetConsoleTextAttribute(hConsoleOut, FOREGROUND_BLUE | FOREGROUND_RED | FOREGROUND_GREEN);
}

///Compares to DNA sequences using the Smith-Waterman method
void CompareDNA::SmithWaterman(DNA *DNASequence1, DNA *DNASequence2,int openGapP, int extendGapP){
	initializeComparisonMatrix(DNASequence1, DNASequence2, 0);
	vector<int> highestCoordinates = scoreSmithWaterman(DNASequence1, DNASequence2, openGapP, extendGapP);
	//printComparisonMatrix();
	traceBackSmithWaterman(DNASequence1, DNASequence2, highestCoordinates);
	calcPercentDifference(DNASequence1, DNASequence2);
}

///Does the scoring for the Smith Waterman method
vector<int> CompareDNA::scoreSmithWaterman(DNA *DNASequence1, DNA *DNASequence2, int openGapP, int extendGapP){
	openingGapPenalty = openGapP;
	vector<int> highestCoordinates(2);
	int highestScore = 0;
	highestCoordinates[1] = 1;
	highestCoordinates[0] = 1;

	int diagonalScore;
	int topScore;
	int leftScore;

	//comparing the two dna sequences and scoring in the vector
	for(int i = 1; i < DNASequence2->getOriginalSize()+1; i++){
		for(int j = 1; j < DNASequence1->getOriginalSize()+1; j++){
			int diagonalI = i - 1;
			int diagonalJ = j - 1;
			int leftI = i;
			int leftJ = j - 1;
			int topI = i - 1;
			int topJ = j;

			//scoring the diagonal top and left locations to the current i,j
			diagonalScore = ComparisonMatrix[diagonalI][diagonalJ] +
					CompareResidue(DNASequence2->getDNAatPositionOriginal(i-1), DNASequence1->getDNAatPositionOriginal(j-1));
			/*if(ComparisonMatrix[topI][topJ] > ComparisonMatrix[leftI][leftJ] &&
			  ComparisonMatrix[topI][topJ] > ComparisonMatrix[diagonalI][diagonalJ]){
				topScore = ComparisonMatrix[topI][topJ] + extendingGapPenalty;
			}else{*/
				topScore = ComparisonMatrix[topI][topJ] + openingGapPenalty;
			//}

			/*if(ComparisonMatrix[leftI][leftJ] > ComparisonMatrix[topI][topJ] &&
			  ComparisonMatrix[leftI][leftJ] > ComparisonMatrix[diagonalI][diagonalJ]){
				leftScore = ComparisonMatrix[leftI][leftJ] + extendingGapPenalty;
			}else{*/
				leftScore = ComparisonMatrix[leftI][leftJ] + openingGapPenalty;
			//}

			//determining the highest score of all and inserting that score into the matrix

			if(0 > diagonalScore && 0 > topScore && 0 > leftScore){
				ComparisonMatrix[i][j] = 0;
			}
			else if(diagonalScore >= topScore && diagonalScore >= leftScore){
				ComparisonMatrix[i][j] = diagonalScore;
			}
			else if(topScore >= diagonalScore && topScore >= leftScore){
				ComparisonMatrix[i][j] = topScore;
			}
			else if(leftScore >= diagonalScore && leftScore >= topScore){
				ComparisonMatrix[i][j] = leftScore;
			}


			//Testing to see if a new high score has been found
			if(ComparisonMatrix[i][j] >= highestScore){
				highestScore = ComparisonMatrix[i][j];
				highestCoordinates[0] = i;
				highestCoordinates[1] = j;
			}
		}
	}
	return highestCoordinates;
}

///Does a trace back of the Matrix created by the Smith Waterman method
void CompareDNA::traceBackSmithWaterman(DNA *DNASequence1, DNA *DNASequence2, vector<int> highestCoordinates){
	//Setting the starting point in the matrix to begin the traceback
	int i = highestCoordinates[0];
	int j = highestCoordinates[1];

	cout << "Begin SmithWaterman traceback" << endl;
	FileIO::Write("Begin SmithWaterman traceback");

	//The actual traceback
	do{
		//creating coordinates for positions around i,j
		int diagonalI = i - 1;
		int diagonalJ = j - 1;
		int leftI = i;
		int leftJ = j - 1;
		int topI = i - 1;
		int topJ = j;

		/*determining the largest neighbor and and depending on which neighbor
		 * it is, setting the aligned sequence of both DNA objects to their respective
		 * residues or gaps. Loop ends when a 0 is the largest neighbor.
		 */
		if(ComparisonMatrix[diagonalI][diagonalJ] > ComparisonMatrix[leftI][leftJ] &&
				ComparisonMatrix[diagonalI][diagonalJ] > ComparisonMatrix[topI][topJ]){
			DNASequence1->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(j-1));
			DNASequence2->InsertAtBeginingAligned(DNASequence2->getDNAatPositionOriginal(i-1));
			i = diagonalI;
			j = diagonalJ;
		}
		else if(ComparisonMatrix[leftI][leftJ] > ComparisonMatrix[topI][topJ] &&
				ComparisonMatrix[leftI][leftJ] > ComparisonMatrix[diagonalI][diagonalJ]){
			DNASequence2->InsertAtBeginingAligned('_');
			DNASequence1->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(j-1));
			j = leftJ;
		}
		else if(ComparisonMatrix[topI][topJ] > ComparisonMatrix[leftI][leftJ] &&
				ComparisonMatrix[topI][topJ] > ComparisonMatrix[diagonalI][diagonalJ]){
				DNASequence2->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(i-1));
				DNASequence1->InsertAtBeginingAligned('_');
				i = topI;
		}
		else{
			DNASequence1->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(j-1));
			DNASequence2->InsertAtBeginingAligned(DNASequence2->getDNAatPositionOriginal(i-1));
			i = 0;
			j = 0;
		}
	}while(ComparisonMatrix[i][j] != 0);
}
/*End Smith-Waterman Functions*/

///Does a comparason using the Needleman/Wunsch method
void CompareDNA::NeedlemanWunsch(DNA *DNASequence1, DNA *DNASequence2, int gapP){
	initializeComparisonMatrix(DNASequence1, DNASequence2, gapP);
	scoreNeedlemanWunsch(DNASequence1, DNASequence2);
	//printScoringMatrix();
	traceBackNeedlemanWunsch(DNASequence1, DNASequence2);
	calcPercentDifference(DNASequence1, DNASequence2);
}

///Responsible for scoring the Needleman/Wunsch method
void CompareDNA::scoreNeedlemanWunsch(DNA *DNASequence1, DNA *DNASequence2){

	//comparing the two dna sequences and scoring in the vector
	for(int i = 1; i < DNASequence2->getOriginalSize()+1; i++){
		for(int j = 1; j < DNASequence1->getOriginalSize()+1; j++){
			int diagonalI = i - 1;
			int diagonalJ = j - 1;
			int leftI = i;
			int leftJ = j - 1;
			int topI = i - 1;
			int topJ = j;

			ComparisonMatrix[i][j] = (ComparisonMatrix[diagonalI][diagonalJ] +
				CompareResidue(DNASequence2->getDNAatPositionOriginal(i-1), DNASequence1->getDNAatPositionOriginal(j-1)));

			if((ComparisonMatrix[leftI][leftJ] + openingGapPenalty) > ComparisonMatrix[i][j]){
				ComparisonMatrix[i][j] = ComparisonMatrix[leftI][leftJ] + openingGapPenalty;
			}
			if((ComparisonMatrix[topI][topJ] + openingGapPenalty) > ComparisonMatrix[i][j]){
				ComparisonMatrix[i][j] = ComparisonMatrix[topI][topJ] + openingGapPenalty;
			}
		}
	}
}

///Preforms a traceback for the matrix created by the Needleman/Wunsch
void CompareDNA::traceBackNeedlemanWunsch(DNA *DNASequence1, DNA *DNASequence2){
	//Finding path and building aligned sequences
	int i = DNASequence2->getOriginalSize();
	int j = DNASequence1->getOriginalSize();
	cout << "Begin actual traceback" << endl;
        FileIO::Write("Begin actual traceback");
	do{
		int diagonalI = i - 1;
		int diagonalJ = j - 1;
		int leftI = i;
		int leftJ = j - 1;
		int topI = i - 1;
		int topJ = j;
		if(DNASequence1->getDNAatPositionOriginal(j-1) == DNASequence2->getDNAatPositionOriginal(i-1)){
			DNASequence1->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(j-1));
			DNASequence2->InsertAtBeginingAligned(DNASequence2->getDNAatPositionOriginal(i-1));
			i = diagonalI;
			j = diagonalJ;
		}
		else if(ComparisonMatrix[topI][topJ] >= ComparisonMatrix[leftI][leftJ]){
			DNASequence1->InsertAtBeginingAligned('_');
			DNASequence2->InsertAtBeginingAligned(DNASequence2->getDNAatPositionOriginal(i-1));
			i = topI;
		}
		else if(ComparisonMatrix[leftI][leftJ] >= ComparisonMatrix[topI][topJ]){
			DNASequence2->InsertAtBeginingAligned('_');
			DNASequence1->InsertAtBeginingAligned(DNASequence1->getDNAatPositionOriginal(j-1));
			j = leftJ;
		}
	}while(i > 0 && j > 0);
}
//End Needleman/Wunsch Functions
